#include "CUDA_HOTFiberTracking.h"

#include "../DWMRICore/DTI.h"
#include "../DWMRICore/HOT.h"
#include "../DWMRICore/Scalar.h"

#include "CUDA_TIJK.h"

/* DEFINES */
#define CUDA_TEN_T2M(m, t) ( \
   (m)[0] = (t)[1], (m)[1] = (t)[2], (m)[2] = (t)[3], \
   (m)[3] = (t)[2], (m)[4] = (t)[4], (m)[5] = (t)[5], \
   (m)[6] = (t)[3], (m)[7] = (t)[5], (m)[8] = (t)[6] )

__shared__ float3 sample_points[SAMPLE_POINT_NUM];
__shared__ float sample_point_weights[SAMPLE_POINT_NUM];

__shared__ float stepSize;
__shared__ float maxLength;
__shared__ float faThreshold;
__shared__ float angleThreshold;
__shared__ float MLSSize;
__shared__ int innerSteps;

/* textures */
/* MASK */
texture<float, cudaTextureType3D> texMASK;
/* DTI */
texture<float4, cudaTextureType3D> texDTI_1;
texture<float4, cudaTextureType3D> texDTI_2;
/* HOT */
texture<float4, cudaTextureType3D> texHOT_1;
texture<float4, cudaTextureType3D> texHOT_2;
texture<float4, cudaTextureType3D> texHOT_3;
texture<float4, cudaTextureType3D> texHOT_4;
texture<float4, cudaTextureType3D> texHOT_5;
texture<float4, cudaTextureType3D> texHOT_6;
texture<float4, cudaTextureType3D> texHOT_7;


/************* TIJK code in CUDA *************/
__device__ const unsigned int _tijk_max_candidates_3d = 30;
__device__ static float _candidates_3d_f[90] = { 
	-0.546405f, 0.619202f, 0.563943f,     
	-0.398931f,-0.600006f, 0.693432f,      
	0.587973f, 0.521686f, 0.618168f,      
	0.055894f,-0.991971f,-0.113444f,       
	-0.666933f,-0.677984f, 0.309094f,      
	0.163684f, 0.533013f, 0.830123f,      
	0.542826f, 0.133898f, 0.829102f,      
	-0.074751f,-0.350412f, 0.933608f,      
	0.845751f, -0.478624f,-0.235847f,      
	0.767148f,-0.610673f, 0.196372f,      
	-0.283810f, 0.381633f, 0.879663f,     
	0.537228f,-0.616249f, 0.575868f,      
	-0.711387f, 0.197778f, 0.674398f,     
	0.886511f, 0.219025f, 0.407586f,      
	0.296061f, 0.842985f, 0.449136f,      
	-0.937540f,-0.340990f, 0.068877f,      
	0.398833f, 0.917023f, 0.000835f,      
	0.097278f,-0.711949f, 0.695460f,      
	-0.311534f, 0.908623f,-0.278121f,      
	-0.432043f,-0.089758f, 0.897375f,      
	-0.949980f, 0.030810f, 0.310788f,     
	0.146722f,-0.811981f,-0.564942f,       
	-0.172201f,-0.908573f, 0.380580f,      
	0.507209f,-0.848578f,-0.150513f,       
	-0.730808f,-0.654136f,-0.194999f,      
	0.077744f, 0.094961f, 0.992441f,      
	0.383976f,-0.293959f, 0.875300f,      
	0.788208f,-0.213656f, 0.577130f,      
	-0.752333f,-0.301447f, 0.585769f,      
	-0.975732f, 0.165497f,-0.143382f      
};

__shared__ unsigned int tijk_max_candidates_3d;
__shared__ float candidates_3d_f[90];

__device__ void cuda_tijk_add_f (float *res, const float *A, const float *B,
	const cuda_tijk_type *type) 
{
	unsigned int i;
	for (i=0; i<type->num; i++)
		*(res++)=*(A++)+*(B++);
}

__device__ void cuda_tijk_sub_f (float *res, const float *A, const float *B,
	const cuda_tijk_type *type) 
{
	unsigned int i;
	for (i=0; i<type->num; i++)
		*(res++)=*(A++)-*(B++);
}

__device__ void cuda_tijk_incr_f (float *res, const float *A, const cuda_tijk_type *type) 
{
	unsigned int i;
	for (i=0; i<type->num; i++)
		*(res++)+=*(A++);
}

__device__ void cuda_tijk_negate_f (float *res, const float *A, const cuda_tijk_type *type) 
{
	unsigned int i;
	for (i=0; i<type->num; i++)
		*(res++)=-*(A++);
}

__device__ void cuda_tijk_scale_f (float *res, const float s, const float *A, const cuda_tijk_type *type) 
{
	unsigned int i;
	for (i=0; i<type->num; i++)
		*(res++)=s*(*A++);
}

__device__ void cuda_tijk_zero_f (float *res, const cuda_tijk_type *type)
{
	unsigned int i;
	for (i=0; i<type->num; i++)
		*(res++)=0.0;
}

__device__ void cuda_tijk_copy_f (float *res, const float *A, const cuda_tijk_type *type) 
{
	unsigned int i;
	for (i=0; i<type->num; i++)
		*(res++)=*(A++);
}
__device__ float _cuda_tijk_6o3d_sym_s_form_f (const float *A, const float *v)
{
	float v00=v[0]*v[0], v01=v[0]*v[1], v02=v[0]*v[2],
		v11=v[1]*v[1], v12=v[1]*v[2], v22=v[2]*v[2];
	return A[0]*v00*v00*v00+
		A[21]*v11*v11*v11+
		A[27]*v22*v22*v22+
		6*(A[1]*v00*v00*v01+
		A[2]*v00*v00*v02+
		A[15]*v01*v11*v11+
		A[20]*v02*v22*v22+
		A[22]*v11*v11*v12+
		A[26]*v12*v22*v22)+
		15*(A[3]*v00*v00*v11+
		A[5]*v00*v00*v22+
		A[10]*v00*v11*v11+
		A[14]*v00*v22*v22+
		A[23]*v11*v11*v22+
		A[25]*v11*v22*v22)+
		30*(A[4]*v00*v00*v12+
		A[16]*v01*v11*v12+
		A[19]*v01*v22*v22)+
		20*(A[6]*v00*v01*v11+
		A[9]*v00*v02*v22+
		A[24]*v11*v12*v22)+
		60*(A[7]*v00*v01*v12+
		A[8]*v00*v01*v22+
		A[11]*v00*v11*v12+
		A[13]*v00*v12*v22+
		A[17]*v01*v11*v22+
		A[18]*v01*v12*v22)+
		90*A[12]*v00*v11*v22;
}

__device__ void _cuda_tijk_6o3d_sym_v_form_f (float *res, const float *A, const float *v)
{
	float v00=v[0]*v[0], v01=v[0]*v[1], v02=v[0]*v[2],                
		v11=v[1]*v[1], v12=v[1]*v[2], v22=v[2]*v[2];                   
	float v00000=v00*v00*v[0], v00001=v00*v00*v[1], v00002=v00*v00*v[2],
		v00011=v00*v01*v[1],  v00012=v00*v01*v[2],  v00022=v00*v02*v[2],  
		v00111=v00*v11*v[1],  v00112=v00*v11*v[2],  v00122=v00*v12*v[2],  
		v00222=v00*v22*v[2],  v01111=v01*v11*v[1],  v01112=v01*v11*v[2],  
		v01122=v01*v12*v[2],  v01222=v01*v22*v[2],  v02222=v02*v22*v[2],  
		v11111=v11*v11*v[1],  v11112=v11*v11*v[2],  v11122=v11*v12*v[2],  
		v11222=v11*v22*v[2],  v12222=v12*v22*v[2],  v22222=v22*v22*v[2];  
	res[0] = A[0]*v00000+                                          
		5*A[1]*v00001+                                              
		5*A[2]*v00002+                                              
		10*A[3]*v00011+                                             
		20*A[4]*v00012+                                             
		10*A[5]*v00022+                                             
		10*A[6]*v00111+                                             
		30*A[7]*v00112+                                             
		30*A[8]*v00122+                                             
		10*A[9]*v00222+                                             
		5*A[10]*v01111+                                             
		20*A[11]*v01112+                                            
		30*A[12]*v01122+                                            
		20*A[13]*v01222+                                            
		5*A[14]*v02222+                                             
		A[15]*v11111+                                               
		5*A[16]*v11112+                                             
		10*A[17]*v11122+                                            
		10*A[18]*v11222+                                            
		5*A[19]*v12222+                                             
		A[20]*v22222;                                               
	res[1] = A[1]*v00000+                                          
		5*A[3]*v00001+                                              
		5*A[4]*v00002+                                              
		10*A[6]*v00011+                                             
		20*A[7]*v00012+                                             
		10*A[8]*v00022+                                             
		10*A[10]*v00111+                                            
		30*A[11]*v00112+                                            
		30*A[12]*v00122+                                            
		10*A[13]*v00222+                                            
		5*A[15]*v01111+                                             
		20*A[16]*v01112+                                            
		30*A[17]*v01122+                                            
		20*A[18]*v01222+                                            
		5*A[19]*v02222+                                             
		A[21]*v11111+                                               
		5*A[22]*v11112+                                             
		10*A[23]*v11122+                                            
		10*A[24]*v11222+                                            
		5*A[25]*v12222+                                             
		A[26]*v22222;                                               
	res[2] = A[2]*v00000+                                          
		5*A[4]*v00001+                                              
		5*A[5]*v00002+                                              
		10*A[7]*v00011+                                             
		20*A[8]*v00012+                                             
		10*A[9]*v00022+                                             
		10*A[11]*v00111+                                            
		30*A[12]*v00112+                                            
		30*A[13]*v00122+                                            
		10*A[14]*v00222+                                            
		5*A[16]*v01111+                                             
		20*A[17]*v01112+                                            
		30*A[18]*v01122+                                            
		20*A[19]*v01222+                                            
		5*A[20]*v02222+                                             
		A[22]*v11111+                                               
		5*A[23]*v11112+                                             
		10*A[24]*v11122+                                            
		10*A[25]*v11222+                                            
		5*A[26]*v12222+                                             
		A[27]*v22222;                                               
}

__device__ float _cuda_tijk_6o3d_sym_mean_f (const float *A) 
{
	return (A[0]+A[21]+A[27]+
		3.0f*(A[3]+A[5]+A[10]+A[14]+A[23]+A[25])+
		6.0f*A[12])/7.0f;
}

__device__ void _cuda_tijk_6o3d_sym_make_iso_f (float *res, const float s) 
{
	/* initialize to zero, then set non-zero elements */
	unsigned int i;
	for (i=0; i<28; i++)
		res[i]=0;
	res[0]=res[21]=res[27]=s;
	res[3]=res[5]=res[10]=res[14]=res[23]=res[25]=0.2f*s;
	res[12]=s/15.0f;
}

__device__ float _cuda_tijk_6o3d_sym_norm_f (const float *A) {
  return sqrtf(_CUDA_TIJK_6O3D_SYM_TSP(A,A));
}

__device__ void _cuda_tijk_6o3d_sym_grad_f (float *res, const float *A, const float *v) {
  float proj, projv[3];
  _cuda_tijk_6o3d_sym_v_form_f (res, A, v);
  CUDA_ELL_3V_SCALE(res,6.0f,res);
  proj=CUDA_ELL_3V_DOT(res,v);
  CUDA_ELL_3V_SCALE(projv,-proj,v);
  CUDA_ELL_3V_INCR(res,projv);
}

__device__ void _cuda_tijk_6o3d_sym_make_rank1_f (float *res, const float s, const float *v) 
{
	float v00=v[0]*v[0], v01=v[0]*v[1], v02=v[0]*v[2],                
		v11=v[1]*v[1], v12=v[1]*v[2], v22=v[2]*v[2];

	res[0]=s*v00*v00*v00; res[1]=s*v00*v00*v01; res[2]=s*v00*v00*v02;   
	res[3]=s*v00*v00*v11; res[4]=s*v00*v00*v12; res[5]=s*v00*v00*v22;   
	res[6]=s*v00*v01*v11; res[7]=s*v00*v01*v12; res[8]=s*v00*v01*v22;   
	res[9]=s*v00*v02*v22; res[10]=s*v00*v11*v11; res[11]=s*v00*v11*v12; 
	res[12]=s*v00*v11*v22; res[13]=s*v00*v12*v22; res[14]=s*v00*v22*v22;
	res[15]=s*v01*v11*v11; res[16]=s*v01*v11*v12; res[17]=s*v01*v11*v22;
	res[18]=s*v01*v12*v22; res[19]=s*v01*v22*v22; res[20]=s*v02*v22*v22;
	res[21]=s*v11*v11*v11; res[22]=s*v11*v11*v12; res[23]=s*v11*v11*v22;
	res[24]=s*v11*v12*v22; res[25]=s*v11*v22*v22; res[26]=s*v12*v22*v22;
	res[27]=s*v22*v22*v22;                                         
}



__device__ int cuda_tijk_init_rank1_3d_f(float *s, float *v, const float *ten, 
	const cuda_tijk_type *type)
{  
	float absmax = -1.0f;                                       
	unsigned int i;                                       
	float *candidate =candidates_3d_f;              
	/*if (type->dim!=3)
		return 1;*/                                          
	for (i=0; i<tijk_max_candidates_3d; i++) {         
		float val = _cuda_tijk_6o3d_sym_s_form_f(ten, candidate);    
		float absval=fabs(val);                               
		if (absval>absmax) {                                
			absmax=absval;                                    
			*s=val;                                          
			CUDA_ELL_3V_COPY(v, candidate);                    
		}                                                 
		candidate += 3;                                     
	}                                                   
	return 0;                                            
}

__device__ int cuda_tijk_init_max_3d_f(float *s, float *v, const float *ten, const cuda_tijk_type *type) 
{
	float max=0;                                            
	unsigned int i;                                        
	float *candidate=candidates_3d_f;
	/*if (type->dim!=3)
		return 1;*/                                            
	*s = max = _cuda_tijk_6o3d_sym_s_form_f(ten, candidate);
	CUDA_ELL_3V_COPY(v, candidate);                          
	for (i=1; i<tijk_max_candidates_3d; i++) {          
		float val;                                            
		candidate += 3;  
		val=_cuda_tijk_6o3d_sym_s_form_f(ten, candidate);           
		if (val>max) {                                       
			max=val;                                           
			*s=val;                                            
			CUDA_ELL_3V_COPY(v, candidate);                      
		}                                                   
	}                                                     
	return 0;                            
}

__device__ int _cuda_tijk_refine_rank1ormax_3d_f(float *s, float *v, const float *ten,
	const cuda_tijk_type *type, const int refinemax) 
{ 
	float isoten[CUDA_TIJK_TYPE_MAX_NUM], anisoten[CUDA_TIJK_TYPE_MAX_NUM];     
	float der[3], iso, anisonorm, anisonorminv, oldval;            
	char sign=(refinemax || *s>0)?1:-1;                           
	float alpha, beta;                                           
	/*if (type->dim!=3)                        
		return 1;*/                                                

	/* It's easier to do the optimization on the deviatoric */       
	iso=_cuda_tijk_6o3d_sym_mean_f(ten);                            
	_cuda_tijk_6o3d_sym_make_iso_f(isoten, iso);                     
	cuda_tijk_sub_f(anisoten, ten, isoten, type);                   
	anisonorm=_cuda_tijk_6o3d_sym_norm_f(anisoten);                       
	if (anisonorm < _CUDA_TIJK_RANK1_eps_start) {                              
		return 0; /* nothing to do */                              
	} else {                                                   
		anisonorminv = 1.0f / anisonorm;                                
	}                                                         
	alpha = beta = _CUDA_TIJK_RANK1_beta * anisonorminv;                           
	oldval = *s - iso;                                              
	/* set initial derivative */                                 
	_cuda_tijk_6o3d_sym_grad_f(der, anisoten, v);                    
	while (1) { /* refine until convergence */                     
		/* stepsize needs to depend on norm to make the descent */     
		/* scale invariant */                                     
		unsigned int armijoct=0;                                   
		float testv[3], val;                                     
		float dist, derlen=CUDA_ELL_3V_LEN(der);                     
		/* determine stepsize based on Armijo's rule */               
		while (1) {                                              
			++armijoct;                                            
			if (armijoct > _CUDA_TIJK_RANK1_maxtry) {                             
				/* failed to find a valid stepsize */                   
				return 2;                                            
			}                                                     
			CUDA_ELL_3V_SCALE_ADD2(testv,1.0f,v,(float)sign*alpha,der);        
			CUDA_ELL_3V_NORM(testv,testv,dist);                      
			dist = 1.0f - CUDA_ELL_3V_DOT(v,testv);                         
			val=_cuda_tijk_6o3d_sym_s_form_f(anisoten,testv);             
			if (sign*val>=sign*oldval+_CUDA_TIJK_RANK1_sigma*derlen*dist) {        
				/* accept step */                                     
				CUDA_ELL_3V_COPY(v,testv);                            
				*s=val+iso;                                          
				_cuda_tijk_6o3d_sym_grad_f(der, anisoten, v);              
				if (alpha<beta) alpha /= _CUDA_TIJK_RANK1_gamma;                   
				break;                                               
			}                                                     
			alpha *= _CUDA_TIJK_RANK1_gamma; /* try again with decreased stepsize */
		}                                                       
		if (sign*(val-oldval)<=_CUDA_TIJK_RANK1_eps_impr*anisonorm) {            
			break; /* declare convergence */                          
		}                                                       
		oldval=val;                                              
	}                                                         
	return 0;                                                  
}

__device__ int cuda_tijk_refine_rank1_3d_f(float *s, float *v, const float *ten, 
	const cuda_tijk_type *type)
{
	return _cuda_tijk_refine_rank1ormax_3d_f(s,v,ten,type,0);
}

__device__ int cuda_tijk_refine_rankk_3d_f(float *ls, float *vs, float *tens,    
	float *res, float *resnorm, const float orignorm,                 
	const cuda_tijk_type *type, const unsigned int k)  
{                                                             
	float newnorm=(*resnorm);                                        
	unsigned int i;                                                 
	/*if (type->dim!=3)                              
		return 1;*/                                                  
	if (*resnorm<_CUDA_TIJK_RANKK_eps_res || k==0) {                         
		return 0; /* nothing to do */                            
	}                                                         
	do {                                                      
		*resnorm=newnorm;                                       
		for (i=0; i<k; i++) {                                   
			if (ls[i]!=0.0) { /* refine an existing term */        
				cuda_tijk_incr_f(res, tens+i*type->num, type);          
				ls[i]=_cuda_tijk_6o3d_sym_s_form_f(res, vs+3*i);         
				if (_CUDA_TIJK_RANKK_pos && ls[i]<0.0f) { /* try a new one */  
					cuda_tijk_init_max_3d_f(ls+i, vs+3*i, res, type);   
				}                                             
			} else { /* add a new term */                        
				if (_CUDA_TIJK_RANKK_pos)                                      
					cuda_tijk_init_max_3d_f(ls+i, vs+3*i, res, type); 
				else                                             
					cuda_tijk_init_rank1_3d_f(ls+i, vs+3*i, res, type);
			}                                                 
			cuda_tijk_refine_rank1_3d_f(ls+i, vs+3*i, res, type);            
			if (!_CUDA_TIJK_RANKK_pos || ls[i]>0.0) {                                
				_cuda_tijk_6o3d_sym_make_rank1_f(tens+i*type->num, ls[i], vs+3*i);
				cuda_tijk_sub_f(res, res, tens+i*type->num, type);  
			} else {                                       
				ls[i]=0.0;                                     
			}                                             
		}                                                 
		newnorm=_cuda_tijk_6o3d_sym_norm_f(res);                        
	} while (newnorm>_CUDA_TIJK_RANKK_eps_res &&                         
		(*resnorm)-newnorm>_CUDA_TIJK_RANKK_eps_impr*orignorm);               
	*resnorm=newnorm;                                       
	return 0;                                              
}


__device__ int cuda_tijk_approx_rankk_3d_f(float *ls, float *vs, float *res, const float *ten,
	const cuda_tijk_type *type, const unsigned int k)  
{                                                                     
/*	char hadls=1, hadvs=1, hadres=1;   */                                                                 
	float newnorm, orignorm;                                     
	unsigned int i;
	float tens[84];
	int retval=0;
	/*if (type->dim!=3)
		return 1;*/                                                
	/* initializations */                                        
	orignorm=newnorm=_cuda_tijk_6o3d_sym_norm_f(ten);                         
	if (orignorm < _CUDA_TIJK_RANKK_eps_res || k == 0) {                                                                  	                                        
		for (i=0; i<k; i++)                                     
			ls[i]=0.0;                                                                                               
		for (int i = 0; i < type->num; ++i) {
			res[i] = ten[i];
		}      
		return 0; /* nothing to do */                               
	}                                                                                                   

	for (i=0; i<k; i++)                                       
		ls[i]=0.0f; 

	for (int i = 0; i < type->num; ++i) {
		res[i] = ten[i];
	}

	retval = cuda_tijk_refine_rankk_3d_f(ls, vs, tens, res, &newnorm,       
		orignorm, type, k);

	return retval;                                              
}

/************* END *************/

__device__ float CUDA_DOT(const float3 *op1, const float3 *op2)
{
	float res = (op1->x) * (op2->x) + (op1->y) * (op2->y) + (op1->z) * (op2->z);
	return res;
}

__device__ float CUDA_GetMASK(const float3 *pos)
{
	float x = pos->x + 0.5f;
	float y = pos->y + 0.5f;
	float z = pos->z + 0.5f;
	float mask = tex3D(texMASK, x, y, z);
	return mask;
}

__device__ void CUDA_GetDTI(const float3 *pos, float *dti)
{
	float x = pos->x + 0.5f;
	float y = pos->y + 0.5f;
	float z = pos->z + 0.5f;

	float4 temp;

	temp = tex3D(texDTI_1, x, y, z);
	dti[0] = temp.x;
	dti[1] = temp.y;
	dti[2] = temp.z;
	dti[3] = temp.w;

	temp = tex3D(texDTI_2, x, y, z);
	dti[4] = temp.x;
	dti[5] = temp.y;
	dti[6] = temp.z;
}

__device__ void CUDA_GetHOT(const float3 *pos, float *hot)
{
	float x = pos->x + 0.5f;
	float y = pos->y + 0.5f;
	float z = pos->z + 0.5f;

	float4 temp;

	temp = tex3D(texHOT_1, x, y, z);
	hot[0] = temp.x;
	hot[1] = temp.y;
	hot[2] = temp.z;
	hot[3] = temp.w;

	temp = tex3D(texHOT_2, x, y, z);
	hot[4] = temp.x;
	hot[5] = temp.y;
	hot[6] = temp.z;
	hot[7] = temp.w;

	temp = tex3D(texHOT_3, x, y, z);
	hot[ 8] = temp.x;
	hot[ 9] = temp.y;
	hot[10] = temp.z;
	hot[11] = temp.w;

	temp = tex3D(texHOT_4, x, y, z);
	hot[12] = temp.x;
	hot[13] = temp.y;
	hot[14] = temp.z;
	hot[15] = temp.w;

	temp = tex3D(texHOT_5, x, y, z);
	hot[16] = temp.x;
	hot[17] = temp.y;
	hot[18] = temp.z;
	hot[19] = temp.w;

	temp = tex3D(texHOT_6, x, y, z);
	hot[20] = temp.x;
	hot[21] = temp.y;
	hot[22] = temp.z;
	hot[23] = temp.w;

	temp = tex3D(texHOT_7, x, y, z);
	hot[24] = temp.x;
	hot[25] = temp.y;
	hot[26] = temp.z;
	hot[27] = temp.w;
}

__device__ float CUDA_DTI2FA(const float *dti)
{
	float t[6];
	t[0] = dti[1]; t[1] = dti[2]; t[2] = dti[3]; t[3] = dti[4]; t[4] = dti[5]; t[5] = dti[6];

	float mean = (t[0] + t[1] + t[2] + t[3] + t[4] + t[5]) / 6.0f;
	if (mean == 0.0f)
		return 0.0f;

	t[0] = t[0] / mean;
	t[1] = t[1] / mean;
	t[2] = t[2] / mean;
	t[3] = t[3] / mean;
	t[4] = t[4] / mean;
	t[5] = t[5] / mean;

	float cross = t[1] * t[1] + t[2] * t[2] + t[4] * t[4];
	float j2 = t[0] * t[3] + t[3] * t[5] + t[5] * t[0] - cross;
	float j4 = t[0] * t[0] + t[3] * t[3] + t[5] * t[5] + 2.0f * cross;
	float fa = sqrtf((j4-j2)/j4);
	return fa;
}

//__device__ int CUDA_DTI2Orientation(const float *dti, float3 *orientation)
//{
//	float t[6];
//	t[0] = dti[1]; t[1] = dti[2]; t[2] = dti[3]; t[3] = dti[4]; t[4] = dti[5]; t[5] = dti[6];
//
//	float mean = (t[0] + t[1] + t[2] + t[3] + t[4] + t[5]) / 6.0f;
//	t[0] = t[0] / mean;
//	t[1] = t[1] / mean;
//	t[2] = t[2] / mean;
//	t[3] = t[3] / mean;
//	t[4] = t[4] / mean;
//	t[5] = t[5] / mean;
//	// comput the major eigen value
//	float I1 = t[0] + t[3] + t[5];
//	float I2 = t[0] * t[3] + t[0] * t[5] + t[3] * t[5] - 
//		(t[1] * t[1] + t[2] * t[2] + t[4] * t[4]);
//	float I3 = t[0] * t[3] * t[5] + 2.0f * (t[1] * t[2] * t[4]) - 
//		(t[0] * t[4] * t[4] + t[3] * t[2] * t[2] + t[5] * t[1] * t[1]);
//
//	float third = 1.0f / 3.0f;
//	float I1third = I1 * third;
//	float I1thirdsqr = I1third * I1third;
//	float I2third = I2 * third;
//	float v = I1thirdsqr - I2third;
//
//	float s = I1thirdsqr * I1third - I1 * I2 / 6.0f + 0.5f * I3;
//
//	if (v < 0.0f)
//		v = 0.0f;
//	float sqrtv = sqrtf(v);
//	float temp = s / (v * sqrtv);
//	if (temp > 1.0f)
//		temp = 1.0f;
//	if (temp < -1.0f)
//		temp = -1.0f;
//	float phi = acosf(temp) * third;
//	float sqrtv2 = 2.0f * sqrtv;
//
//	float e1 = I1third + sqrtv2 * cosf(phi);
//
//	// compute the major eigen vector
//	float3 v1;
//	v1.x = (t[1] * t[4] - (t[3] - e1) * t[2]) * (t[2] * t[4] - (t[5] - e1) * t[1]);
//	v1.y = (t[4] * t[2] - (t[5] - e1) * t[1]) * (t[1] * t[2] - (t[0] - e1) * t[4]);
//	v1.z = (t[4] * t[1] - (t[3] - e1) * t[2]) * (t[2] * t[1] - (t[0] - e1) * t[4]);
//
//	// normalize
//	float len = sqrtf(v1.x * v1.x + v1.y * v1.y + v1.z * v1.z);
//	if (len > 0.0f) {
//		v1.x = v1.x / len;
//		v1.y = v1.y / len;
//		v1.z = v1.z / len;
//	}
//
//	*orientation = v1;
//
//	return 1;
//}

__device__ int CUDA_HOT2Orientations(const float *hot, float3 *orientations)
{
    float ten[28];
    for (int i = 0; i < 28; ++i) {
		ten[i] = hot[i];
	}

	cuda_tijk_type type = {6, 3, 28};

    float ls[4][3], vs[4][9], res[28];
	float p_ten[28];
	float res_norm[3];

	/* first get the positive approximation of the tensor */
	cuda_tijk_approx_rankk_3d_f(ls[3], vs[3], res, ten, &type, 6);
	cuda_tijk_sub_f(p_ten, ten, res, &type);

	float last_norm = _cuda_tijk_6o3d_sym_norm_f(p_ten);

    /* try rank-1 decomposition */
    cuda_tijk_approx_rankk_3d_f(ls[0], vs[0], res, p_ten, &type, 1);
    res_norm[0] = _cuda_tijk_6o3d_sym_norm_f(res);
    if (res_norm[0] < last_norm * HOT_TD_NORM_THRESHOLD) {
		orientations[0].x = vs[0][0];
		orientations[0].y = vs[0][1];
		orientations[0].z = vs[0][2];
        return 1;
    }

    /* then try the rank 2 decomposition */
    cuda_tijk_approx_rankk_3d_f(ls[1], vs[1], res, p_ten, &type, 2);
    res_norm[1] = _cuda_tijk_6o3d_sym_norm_f(res);
    if (res_norm[1] < last_norm * HOT_TD_NORM_THRESHOLD) {
		/* success */
		if (ls[1][1] < ls[1][0] * HOT_TD_LS_THRESHOLD) {
			orientations[0].x = vs[0][0];
			orientations[0].y = vs[0][1];
			orientations[0].z = vs[0][2];
			return 1;
		} else {
			orientations[0].x = vs[1][0];
			orientations[0].y = vs[1][1];
			orientations[0].z = vs[1][2];
			orientations[1].x = vs[1][3];
			orientations[1].y = vs[1][4];
			orientations[1].z = vs[1][5];
			return 2;
		}
    }

    /* finally, use the rank 3 decomposition */
    cuda_tijk_approx_rankk_3d_f(ls[2], vs[2], res, p_ten, &type, 3);
	if (ls[2][1] < ls[2][0] * HOT_TD_LS_THRESHOLD) {
		orientations[0].x = vs[0][0];
		orientations[0].y = vs[0][1];
		orientations[0].z = vs[0][2];
		return 1;
	} else if (ls[2][2] < ls[2][0] * HOT_TD_LS_THRESHOLD) {
		orientations[0].x = vs[1][0];
		orientations[0].y = vs[1][1];
		orientations[0].z = vs[1][2];
		orientations[1].x = vs[1][3];
		orientations[1].y = vs[1][4];
		orientations[1].z = vs[1][5];
		return 2;
	} else {
		orientations[0].x = vs[2][0];
		orientations[0].y = vs[2][1];
		orientations[0].z = vs[2][2];
		orientations[1].x = vs[2][3];
		orientations[1].y = vs[2][4];
		orientations[1].z = vs[2][5];
		orientations[2].x = vs[2][6];
		orientations[2].y = vs[2][7];
		orientations[2].z = vs[2][8];
		return 3;
	}
}



__device__ void CUDA_DTI_Filter(const float3 *pos, float *dti)
{
	float temp_dti[7];
	CUDA_GetDTI(pos, temp_dti);

	float m[9];
	CUDA_TEN_T2M(m, temp_dti);

	/* get the maximum length */
	float max_length = 0.0f;
	for (int i = 0; i < SAMPLE_POINT_NUM; ++i) {
		float3 sample;
		sample.x = sample_points[i].x * m[0] + sample_points[i].y * m[1] + sample_points[i].z * m[2];
		sample.y = sample_points[i].x * m[3] + sample_points[i].y * m[4] + sample_points[i].z * m[5];
		sample.z = sample_points[i].x * m[6] + sample_points[i].y * m[7] + sample_points[i].z * m[8];
		float length = sqrtf(CUDA_DOT(&sample, &sample));
		if (length > max_length)
			max_length = length;
	}

	/* filter */
	for (int i = 0; i < 7; ++i)
		dti[i] = 0.0f;
	float weight = 0.0f;
	for (int i = 0; i < SAMPLE_POINT_NUM; ++i) {
		float3 sample;
		sample.x = sample_points[i].x * m[0] + sample_points[i].y * m[1] + sample_points[i].z * m[2];
		sample.y = sample_points[i].x * m[3] + sample_points[i].y * m[4] + sample_points[i].z * m[5];
		sample.z = sample_points[i].x * m[6] + sample_points[i].y * m[7] + sample_points[i].z * m[8];
		sample.x = sample.x / max_length * MLSSize + pos->x;
		sample.y = sample.y / max_length * MLSSize + pos->y;
		sample.z = sample.z / max_length * MLSSize + pos->z;
		if (CUDA_GetMASK(&sample) < 0.01f) {
			continue;
		}
		CUDA_GetDTI(&sample, temp_dti);
		for (int j = 0; j < 7; ++j) {
			dti[j] += sample_point_weights[i] * temp_dti[j];
		}
		weight += sample_point_weights[i];
	}
	
	for (int j = 0; j < 7; ++j) {
		dti[j] = dti[j] / weight;
	}
}

__device__ void CUDA_HOT_Filter(const float3 *pos, float *hot)
{
	float temp_dti[7];
	CUDA_GetDTI(pos, temp_dti);

	float m[9];
	CUDA_TEN_T2M(m, temp_dti);

	/* get the maximum length */
	float max_length = 0.0f;
	for (int i = 0; i < SAMPLE_POINT_NUM; ++i) {
		float3 sample;
		sample.x = sample_points[i].x * m[0] + sample_points[i].y * m[1] + sample_points[i].z * m[2];
		sample.y = sample_points[i].x * m[3] + sample_points[i].y * m[4] + sample_points[i].z * m[5];
		sample.z = sample_points[i].x * m[6] + sample_points[i].y * m[7] + sample_points[i].z * m[8];
		float length = sqrtf(CUDA_DOT(&sample, &sample));
		if (length > max_length)
			max_length = length;
	}

	float temp_hot[28];
	/* filter */
	for (int i = 0; i < 28; ++i)
		hot[i] = 0.0f;
	float weight = 0.0f;
	for (int i = 0; i < SAMPLE_POINT_NUM; ++i) {
		float3 sample;
		sample.x = sample_points[i].x * m[0] + sample_points[i].y * m[1] + sample_points[i].z * m[2];
		sample.y = sample_points[i].x * m[3] + sample_points[i].y * m[4] + sample_points[i].z * m[5];
		sample.z = sample_points[i].x * m[6] + sample_points[i].y * m[7] + sample_points[i].z * m[8];
		sample.x = sample.x / max_length * MLSSize + pos->x;
		sample.y = sample.y / max_length * MLSSize + pos->y;
		sample.z = sample.z / max_length * MLSSize + pos->z;
		if (CUDA_GetMASK(&sample) < 0.01f) {
			continue;
		}
		CUDA_GetHOT(&sample, temp_hot);
		for (int j = 0; j < 28; ++j) {
			hot[j] += sample_point_weights[i] * temp_hot[j];
		}
		weight += sample_point_weights[i];
	}
	
	for (int j = 0; j < 28; ++j) {
		hot[j] = hot[j] / weight;
	}
}

__device__ int CUDA_GetOrientation(const float3 *in_pos, const float3 *in_dir, float3 *out_dirs)
{
	out_dirs[0].x = out_dirs[0].y = out_dirs[0].z = 0.0f;

	float mask = CUDA_GetMASK(in_pos);
	if (mask < 0.01f)
		return 0;

	float dti[7];
	CUDA_DTI_Filter(in_pos, dti);

	float fa = CUDA_DTI2FA(dti);
	if (fa < faThreshold)
		return 0;

	float hot[28];
	float3 dirs[3];
	int num = CUDA_HOT2Orientations(hot, dirs);

	if (in_dir == NULL) {                       /* no input orientations */
		for (int i = 0; i < num; ++i)
			out_dirs[i] = dirs[i];
		return num;
	} else {                                    /* select the best one */
		float max_v = 0.0f;
		int index = 0;
		for (int i = 0; i < num; ++i) {
			float v = CUDA_DOT(in_dir, &(dirs[i]));
			if (fabs(v) > fabs(max_v)) {
				max_v = v;
				index = i;
			}
		}
		out_dirs[0].x = dirs[index].x;
		out_dirs[0].y = dirs[index].y;
		out_dirs[0].z = dirs[index].z;

		/* check the angle threshold */
		if (fabs(max_v) < angleThreshold)
			return 0;

		/* re-orientation */
		if (max_v < 0.0f) {
			out_dirs[0].x = -(out_dirs[0].x);
			out_dirs[0].y = -(out_dirs[0].y);
			out_dirs[0].z = -(out_dirs[0].z);
		}
		return 1;
	}
}

__device__ void CUDA_UpdatePosition(const float3 *last_pos, const float3 *last_dir, float3 *next_pos)
{
	next_pos->x = last_pos->x + last_dir->x * stepSize;
	next_pos->y = last_pos->y + last_dir->y * stepSize;
	next_pos->z = last_pos->z + last_dir->z * stepSize;
}

__device__ int CUDA_ComputeNextPositionRK2(const float3 *last_pos, const float3 *last_dir,
										   float3 *next_pos, float3* next_dir)
{
	float3 p;
	p.x = last_pos->x + last_dir->x * stepSize / 2.0f;
	p.y = last_pos->y + last_dir->y * stepSize / 2.0f;
	p.z = last_pos->z + last_dir->z * stepSize / 2.0f;

	float3 v;
	if (CUDA_GetOrientation(&p, last_dir, &v) == 0) {
		return 0;
	}

	next_pos->x = last_pos->x + v.x * stepSize;
	next_pos->y = last_pos->y + v.y * stepSize;
	next_pos->z = last_pos->z + v.z * stepSize;

	if (CUDA_GetOrientation(next_pos, last_dir, next_dir) == 0) {
		return 0;
	}

	return 1;
}


__global__ void 
GenerateSeeds_kernel(const float3 *points, const float *weights, CUDA_Seed *seeds, 
					 const float in_faThreshold, const float in_MLSSize, const int count)
{
	int x = blockIdx.x * blockDim.x + threadIdx.x;
	int y = blockIdx.y * blockDim.y + threadIdx.y;

	int index = y * gridDim.x * blockDim.x + x;

	if (threadIdx.x == 0) {
		tijk_max_candidates_3d = _tijk_max_candidates_3d;
		for (int i = 0; i < _tijk_max_candidates_3d; ++i) {
			candidates_3d_f[i*3+0] = _candidates_3d_f[i*3+0];
			candidates_3d_f[i*3+1] = _candidates_3d_f[i*3+1];
			candidates_3d_f[i*3+2] = _candidates_3d_f[i*3+2];
		}
		for (int i = 0; i < SAMPLE_POINT_NUM; ++i) {
			sample_points[i] = points[i];
			sample_point_weights[i] = weights[i];
		}

		 faThreshold = in_faThreshold;
		 MLSSize = in_MLSSize;
	}
	__syncthreads();

	if (index < count) {
		float3 pos = seeds[index].pos;

		float3 dir;
		int num = CUDA_GetOrientation(&pos, NULL, &dir);

		seeds[index].dir = dir;

		if (num == 0)
			seeds[index].id = -1;
	}

	//__syncthreads();
}



__global__ void
TrackFibers_kernel(const float3 *points, const float *weights, CUDA_Seed *seeds, 
				   CUDA_Fiber *fibers, const float in_stepSize, const float in_maxLength, 
				   const float in_faThreshold, const float in_angleThreshold, 
				   const float in_MLSSize, const int in_innerSteps, const int count)
{
	int x = blockIdx.x * blockDim.x + threadIdx.x;
	int y = blockIdx.y * blockDim.y + threadIdx.y;

	int index = y * gridDim.x * blockDim.x + x;

	if (threadIdx.x == 0) {
		tijk_max_candidates_3d = _tijk_max_candidates_3d;
		for (int i = 0; i < _tijk_max_candidates_3d; ++i) {
			candidates_3d_f[i*3+0] = _candidates_3d_f[i*3+0];
			candidates_3d_f[i*3+1] = _candidates_3d_f[i*3+1];
			candidates_3d_f[i*3+2] = _candidates_3d_f[i*3+2];
		}

		for (int i = 0; i < SAMPLE_POINT_NUM; ++i) {
			sample_points[i] = points[i];
			sample_point_weights[i] = weights[i];
		}

		stepSize = in_stepSize;
		maxLength = in_maxLength;
		faThreshold = in_faThreshold;
		angleThreshold = in_angleThreshold;
		innerSteps = in_innerSteps;
		MLSSize = in_MLSSize;
	}
	__syncthreads();

	if (index < count) {
		const CUDA_Seed *seed = &(seeds[index]);
		CUDA_Fiber *fiber = &(fibers[index]);

		int maxSteps = (int)(maxLength / (stepSize * (float)innerSteps));

		int fnum = 0;
		int bnum = 0;
		float3 last_dir = seed->dir;

		/* tracking forward */
		fiber->f[0] = seed->pos;
		for (fnum = 1; fnum < maxSteps; ++fnum) {
			float3 last_pos = fiber->f[fnum-1];
			bool flag = false;
			for (int step = 0; step < innerSteps; ++step) {
				float3 next_pos, next_dir;
				if (CUDA_ComputeNextPositionRK2(&last_pos, &last_dir, &next_pos, &next_dir) == 0) {
					flag = true;
					break;
				}
				last_dir = next_dir;
				last_pos = next_pos;
			}
			if (flag == true)
				break;
			fiber->f[fnum] = last_pos;
		}

		///* tracking backward */
		last_dir.x = -seed->dir.x; last_dir.y = -seed->dir.y; last_dir.z = -seed->dir.z;
		fiber->b[0] = seed->pos;
		for (bnum = 1; bnum < maxSteps; ++bnum) {
			float3 last_pos = fiber->b[bnum-1];
			bool flag = false;
			for (int step = 0; step < innerSteps; ++step) {
				float3 next_pos, next_dir;
				if (CUDA_ComputeNextPositionRK2(&last_pos, &last_dir, &next_pos, &next_dir) == 0) {
					flag = true;
					break;
				}
				last_dir = next_dir;
				last_pos = next_pos;
			}
			if (flag == true)
				break;
			fiber->b[bnum] = last_pos;
		}

		fiber->fc = fnum;
		fiber->bc = bnum;

		fiber->seed.dir = seed->dir;
		fiber->seed.pos = seed->pos;
		fiber->seed.id = seed->id;
	}
}


/************************************************************************/
/************************************************************************/
/************************************************************************/
extern "C"
void CUDA_GenerateSeeds(GPUPlan *plan, const float in_faThreshold, 
						const float in_MLSSize, const int num)
{
	cutilSafeCall(cudaSetDevice(plan->deviceID));
	cutilSafeCall(cudaMemcpyAsync(plan->d_pSeeds, plan->h_pSeeds, 
		THREADS_PER_RUN*sizeof(CUDA_Seed), cudaMemcpyHostToDevice, plan->stream));

	dim3 dimBlock(64, 1);
	dim3 dimGrid(THREADS_PER_RUN/64, 1);
	GenerateSeeds_kernel<<<dimGrid, dimBlock, 0, plan->stream>>>
		(plan->d_pSamplePoints, plan->d_pSampleWeights, plan->d_pSeeds, in_faThreshold, in_MLSSize, num);

	cutilSafeCall(cudaGetLastError());

	cutilSafeCall(cudaMemcpyAsync(plan->h_pSeeds, plan->d_pSeeds, 
		THREADS_PER_RUN*sizeof(CUDA_Seed), cudaMemcpyDeviceToHost, plan->stream));
}

extern "C"
void CUDA_TrackFibers(GPUPlan *plan, const float in_stepSize, const float in_maxLength, 
					  const float in_faThreshold, const float in_angleThreshold, 
					  const float in_MLSSize, const int in_innerSteps, const int num)
{
	cutilSafeCall(cudaSetDevice(plan->deviceID));
	cutilSafeCall(cudaMemcpyAsync(plan->d_pSeeds, plan->h_pSeeds, 
		THREADS_PER_RUN*sizeof(CUDA_Seed), cudaMemcpyHostToDevice, plan->stream));

	dim3 dimBlock(512, 1);
	dim3 dimGrid(THREADS_PER_RUN/512, 1);
	TrackFibers_kernel<<<dimGrid, dimBlock, 0, plan->stream>>>
		(plan->d_pSamplePoints, plan->d_pSampleWeights, plan->d_pSeeds, 
		 plan->d_pFibers, in_stepSize, in_maxLength, in_faThreshold, 
		 in_angleThreshold, in_MLSSize, in_innerSteps, num);


	cutilSafeCall(cudaGetLastError());

	cutilSafeCall(cudaMemcpyAsync(plan->h_pFibers, plan->d_pFibers, 
		THREADS_PER_RUN*sizeof(CUDA_Fiber), cudaMemcpyDeviceToHost, plan->stream));
}

extern "C"
void CUDA_BindTextures(CHOT *HOT, CDTI *DTI, CScalar *MASK, GPUPlan *plans)
{
	int w = DTI->GetWidth();
	int h = DTI->GetHeight();
	int d = DTI->GetDepth();
	int dti_n = DTI->GetN();
	int hot_n = HOT->GetN();

	float *hot = HOT->GetHOTData();
	float *dti = DTI->GetDTIData();
	float *mask = MASK->GetScalarData();

	/* create textures and bind them */

	/* MASK */
	cudaExtent extent = make_cudaExtent(w, h, d);
	cudaChannelFormatDesc desc = cudaCreateChannelDesc<float>();
	cutilSafeCall(cudaMalloc3DArray(&(plans->d_pMASK), &desc, extent));

	cudaMemcpy3DParms copyParams = {0};
	copyParams.srcPtr   = make_cudaPitchedPtr((void*)(mask), w*sizeof(float), w, h);
	copyParams.dstArray = plans->d_pMASK;
	copyParams.extent   = extent;
	copyParams.kind     = cudaMemcpyHostToDevice;
	cutilSafeCall(cudaMemcpy3D(&copyParams));

	texMASK.normalized = false;                    
	texMASK.filterMode = cudaFilterModeLinear;      
	texMASK.addressMode[0] = cudaAddressModeClamp;
	texMASK.addressMode[1] = cudaAddressModeClamp;
	texMASK.addressMode[2] = cudaAddressModeClamp;
	cutilSafeCall(cudaBindTextureToArray(texMASK, plans->d_pMASK, desc));


	/* DTI */
	for (int k = 0; k < 2; ++k){
		desc = cudaCreateChannelDesc<float4>();
		cutilSafeCall(cudaMalloc3DArray(&(plans->d_pDTI[k]), &desc, extent));

		float *p = new float[w*h*d*4];
		for (int j = 0; j < w*h*d; ++j) {
			p[j*4+0] = dti[j*dti_n+k*4+0];
			p[j*4+1] = dti[j*dti_n+k*4+1];
			p[j*4+2] = dti[j*dti_n+k*4+2];
			if (k == 0)
				p[j*4+3] = dti[j*dti_n+k*4+3];
			else
				p[j*4+3] = 0.0f;
		}

		cudaMemcpy3DParms copyParams = {0};
		copyParams.srcPtr   = make_cudaPitchedPtr((void*)(p), w*sizeof(float4), w, h);
		copyParams.dstArray = plans->d_pDTI[k];
		copyParams.extent   = extent;
		copyParams.kind     = cudaMemcpyHostToDevice;
		cutilSafeCall(cudaMemcpy3D(&copyParams));

		if (k == 0) {
			texDTI_1.normalized = false;
			texDTI_1.filterMode = cudaFilterModeLinear;      
			texDTI_1.addressMode[0] = cudaAddressModeClamp;
			texDTI_1.addressMode[1] = cudaAddressModeClamp;
			texDTI_1.addressMode[2] = cudaAddressModeClamp;
			cutilSafeCall(cudaBindTextureToArray(texDTI_1, plans->d_pDTI[k], desc));
		} else {
			texDTI_2.normalized = false;
			texDTI_2.filterMode = cudaFilterModeLinear;      
			texDTI_2.addressMode[0] = cudaAddressModeClamp;
			texDTI_2.addressMode[1] = cudaAddressModeClamp;
			texDTI_2.addressMode[2] = cudaAddressModeClamp;
			cutilSafeCall(cudaBindTextureToArray(texDTI_2, plans->d_pDTI[k], desc));
		}

		delete[] p;
	}

	/* HOT */
	for (int k = 0; k < 7; ++k){
		desc = cudaCreateChannelDesc<float4>();
		cutilSafeCall(cudaMalloc3DArray(&(plans->d_pHOT[k]), &desc, extent));

		float *p = new float[w*h*d*4];
		for (int j = 0; j < w*h*d; ++j) {
			p[j*4+0] = hot[j*hot_n+k*4+0];
			p[j*4+1] = hot[j*hot_n+k*4+1];
			p[j*4+2] = hot[j*hot_n+k*4+2];
			p[j*4+3] = hot[j*hot_n+k*4+3];
		}

		cudaMemcpy3DParms copyParams = {0};
		copyParams.srcPtr   = make_cudaPitchedPtr((void*)(p), w*sizeof(float4), w, h);
		copyParams.dstArray = plans->d_pHOT[k];
		copyParams.extent   = extent;
		copyParams.kind     = cudaMemcpyHostToDevice;
		cutilSafeCall(cudaMemcpy3D(&copyParams));

		if (k == 0) {
			texHOT_1.normalized = false;
			texHOT_1.filterMode = cudaFilterModeLinear;      
			texHOT_1.addressMode[0] = cudaAddressModeClamp;
			texHOT_1.addressMode[1] = cudaAddressModeClamp;
			texHOT_1.addressMode[2] = cudaAddressModeClamp;
			cutilSafeCall(cudaBindTextureToArray(texHOT_1, plans->d_pHOT[k], desc));
		} else if(k == 1) {
			texHOT_2.normalized = false;
			texHOT_2.filterMode = cudaFilterModeLinear;      
			texHOT_2.addressMode[0] = cudaAddressModeClamp;
			texHOT_2.addressMode[1] = cudaAddressModeClamp;
			texHOT_2.addressMode[2] = cudaAddressModeClamp;
			cutilSafeCall(cudaBindTextureToArray(texHOT_2, plans->d_pHOT[k], desc));
		} else if (k == 2) {
			texHOT_3.normalized = false;
			texHOT_3.filterMode = cudaFilterModeLinear;      
			texHOT_3.addressMode[0] = cudaAddressModeClamp;
			texHOT_3.addressMode[1] = cudaAddressModeClamp;
			texHOT_3.addressMode[2] = cudaAddressModeClamp;
			cutilSafeCall(cudaBindTextureToArray(texHOT_3, plans->d_pHOT[k], desc));
		} else if (k == 3) {
			texHOT_4.normalized = false;
			texHOT_4.filterMode = cudaFilterModeLinear;      
			texHOT_4.addressMode[0] = cudaAddressModeClamp;
			texHOT_4.addressMode[1] = cudaAddressModeClamp;
			texHOT_4.addressMode[2] = cudaAddressModeClamp;
			cutilSafeCall(cudaBindTextureToArray(texHOT_4, plans->d_pHOT[k], desc));
		} else if (k == 4) {
			texHOT_5.normalized = false;
			texHOT_5.filterMode = cudaFilterModeLinear;      
			texHOT_5.addressMode[0] = cudaAddressModeClamp;
			texHOT_5.addressMode[1] = cudaAddressModeClamp;
			texHOT_5.addressMode[2] = cudaAddressModeClamp;
			cutilSafeCall(cudaBindTextureToArray(texHOT_5, plans->d_pHOT[k], desc));
		} else if (k == 5) {
			texHOT_6.normalized = false;
			texHOT_6.filterMode = cudaFilterModeLinear;      
			texHOT_6.addressMode[0] = cudaAddressModeClamp;
			texHOT_6.addressMode[1] = cudaAddressModeClamp;
			texHOT_6.addressMode[2] = cudaAddressModeClamp;
			cutilSafeCall(cudaBindTextureToArray(texHOT_6, plans->d_pHOT[k], desc));
		} else {
			texHOT_7.normalized = false;
			texHOT_7.filterMode = cudaFilterModeLinear;      
			texHOT_7.addressMode[0] = cudaAddressModeClamp;
			texHOT_7.addressMode[1] = cudaAddressModeClamp;
			texHOT_7.addressMode[2] = cudaAddressModeClamp;
			cutilSafeCall(cudaBindTextureToArray(texHOT_7, plans->d_pHOT[k], desc));
		}

		delete[] p;
	}
}
